In [1]:

    

%pylab inline
figsize(15, 12)

import os
from functools import partial

from skimage import io, measure, morphology
from scipy import ndimage as nd
import mahotas as mh


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

cd "/Users/nuneziglesiasj/Dropbox/RNAPII Quant Juan/Telomere tri-colour staining"


    

    




/Users/nuneziglesiasj/Dropbox/rnapii quant juan/Telomere tri-colour staining

In [3]:

    

samples = filter(lambda x: x.endswith('tif_files'), os.listdir('.'))
images = []
for s in samples:
    im_fns = filter(lambda x: x.endswith('.tif'), os.listdir(s))
    im_fns = [os.path.join(s, im_fn) for im_fn in im_fns]
    ims = map(mh.imread, im_fns)
    ims = [im[..., 0] for im in ims] # images are read as grayscale 3-channel RGB. Keep only first channel
    images.append(ims)


    

In [4]:

    

centromeres = [chs[1] for chs in images]


    

In [5]:

    

from matplotlib import cm
def imshowa(im):
    imshow(im, cmap=cm.cubehelix, interpolation='nearest')


    

In [6]:

    

def imshow_grid(imlist):
    for i in range(5):
        subplot(231+i)
        imshowa(imlist[i])
        xticks([]); yticks([])
    tight_layout()

imshow_grid(centromeres)


    

In [7]:

    

from skimage import filter as imfilter
centromere_ts_otsu = map(imfilter.threshold_otsu, centromeres)
centromere_thresholded_otsu = map(lambda x: x[0] > x[1], zip(centromeres, centromere_ts_otsu))


    

In [8]:

    

imshow_grid(centromere_thresholded_otsu)


    

In [9]:

    

centromere_thresholded_noiseless = map(partial(morphology.remove_small_objects, min_size=36), centromere_thresholded_otsu)
centromere_radius = 10
centromere_neighborhood = map(partial(nd.binary_dilation, iterations=centromere_radius), centromere_thresholded_noiseless)
imshow_grid(centromere_neighborhood)


    

In [10]:

    

telomeres = [chs[2] for chs in images]


    

In [11]:

    

imshow_grid(telomeres)


    

In [12]:

    

telomeres_threshold_adapt = map(partial(imfilter.threshold_adaptive, block_size=49), telomeres)
telomeres_threshold_adapt = map(partial(nd.binary_opening, structure=morphology.disk(4)), telomeres_threshold_adapt)
imshow_grid(telomeres_threshold_adapt)


    

In [13]:

    

subplot(121)
imshow(telomeres[0], cmap=cm.cubehelix, interpolation='nearest')
subplot(122)
imshow(telomeres_threshold_adapt[0], interpolation='nearest')
tight_layout()


    

In [14]:

    

non_centromere_telomeres = map( # lambda x: x[1] - (x[0] * x[1])
                               lambda x: 2 * x[0].astype(uint8) + x[1], zip(centromere_thresholded_noiseless, telomeres_threshold_adapt))
imshow_grid(non_centromere_telomeres)


    

In [15]:

    

import cafe
encoded0 = cafe.encode_centro_telomeres(centromeres[0], telomeres[0])
imshowa(encoded0)


    

In [16]:

    

from skimage import data
from skimage import viewer
from skimage.viewer.plugins.overlayplugin import OverlayPlugin
from skimage.viewer.widgets import Slider
# how to use the viewer:
# >>> v = viewer.ImageViewer(data.camera())
# >>> v.show()


    

In [23]:

    

def encode(im, **kwargs):
    if kwargs.has_key('alpha'):
        a = kwargs['alpha']
        del kwargs['alpha']
    else:
        a = 100
    ov = cafe.encode_centro_telomeres_multichannel(im, **kwargs)
    return (a * (ov > 0) + int(a * 0.42) * ov).astype(np.uint8)

class CentroPlugin(OverlayPlugin):
    def __init__(self, *args, **kwargs):
        super(CentroPlugin, self).__init__(image_filter=encode, **kwargs)
        
    def attach(self, image_viewer):
        self.add_widget(Slider('alpha', 0, 100, value=100, value_type='int'))
        self.add_widget(Slider('centro_min_size', 0, 100, value=14, value_type='int'))
        self.add_widget(Slider('centro_radius', 0, 100, value=10, value_type='int'))
        self.add_widget(Slider('telo_offset', -256, 255, value=0))
        self.add_widget(Slider('telo_adapt_radius', 0, 101, value=49, value_type='int'))
        self.add_widget(Slider('telo_open_radius', 0, 20, value=4, value_type='int'))
        super(CentroPlugin, self).attach(image_viewer)


    

In [18]:

    

dapis = [chs[0] for chs in images]
rgbs = map(np.dstack, zip(centromeres, telomeres, dapis))


    

In [19]:

    

reload(cafe)


    

    
Out[19]:





<module 'skimage.viewer' from '/Users/nuneziglesiasj/venv/cafe/lib/python2.7/site-packages/scikit_image-0.10dev-py2.7-macosx-10.5-x86_64.egg/skimage/viewer/__init__.pyc'>

In [20]:

    

v = viewer.ImageViewer(rgbs[0])
v += CentroPlugin()
overlay = v.show()[0][0]


    

In [27]:

    

from skimage import segmentation as seg
overlay = seg.relabel_sequential(overlay)[0]
imshowa(overlay)


    

In [30]:

    

aukbs = [chs[3] for chs in images]
masks = []
measurements = []
for s, aukb, rgb in zip(samples, aukbs, rgbs):
    v = viewer.ImageViewer(rgb)
    v += CentroPlugin()
    overlay = v.show()[0][0]
    overlay = seg.relabel_sequential(overlay)[0]
    mask = (overlay == 1)
    masks.append(mask)
    aukb_measurement = aukb[mask]
    measurements.append(aukb_measurement)
    fout_txt = os.path.join(s, 'measure.txt')
    np.savetxt(fout_txt, aukb_measurement)
    fout_im = os.path.join(s, 'mask.png')
    mh.imsave(fout_im, 64 * overlay.astype(np.uint8))


    

In [35]:

    

plt.hist(measurements[4], bins=25)


    

    
Out[35]:





(array([  38.,   80.,  193.,  283.,  498.,  373.,  291.,  252.,  217.,
        166.,  148.,   90.,   74.,   56.,   45.,   33.,   26.,   21.,
         11.,   22.,    9.,    8.,    4.,    4.,    3.]),
 array([   4.  ,    7.96,   11.92,   15.88,   19.84,   23.8 ,   27.76,
         31.72,   35.68,   39.64,   43.6 ,   47.56,   51.52,   55.48,
         59.44,   63.4 ,   67.36,   71.32,   75.28,   79.24,   83.2 ,
         87.16,   91.12,   95.08,   99.04,  103.  ]),
 <a list of 25 Patch objects>)